rm(list=ls())

require(cjoint)
require(tidyverse)
library(extrafont)
library(ashr)

# load
dat <- read_csv("data/cj_data.csv")
dat <- na.exclude(dat)

str(dat)
dat[,c(43:45)] <- data.frame(lapply(dat[,c(43:45)],as.factor))
dat[,c(47:48)] <- data.frame(lapply(dat[,c(47:48)],as.factor))

dat <- dat |>
  dplyr::select(- ...1, respondent)

# AMCE whole
attr_list <- list()             
attr_list[["General_description"]] <- c("De fact legalization","Right to choose")
attr_list[["Religious_customary"]] <- c("Religious laws", "Tribal customs")
attr_list[["Support_criticism"]] <- c("Legal fragmentaion", "Protect poor women", "Social order", "Violation to women's right")
attr_list[["Internal_external"]] <- c("Domestic criticism", "International criticism", "Evaluation")

cjoint_design <- makeDesign(type = "constraints", 
                            attribute.levels = attr_list)

cjoint_pool <- amce(selected ~ General_description + Religious_customary + Support_criticism + Internal_external, 
                    data = dat)

plot(cjoint_pool, point.size = 1.2, dodge.size = 1.2, text.size = 18)
summary(cjoint_pool)


# Bonferroni correction
res <- summary(cjoint_pool)$amce

point.e <- res$Estimate
sd.errors <- res$`Std. Err`

names(point.e) <- paste(res$Attribute, res$Level)
names(sd.errors) <- paste(res$Attribute, res$Level)

p.raw <- 2 * (1 - pnorm(abs(point.e / sd.errors)))

p.bonf <- p.adjust(p.raw, method = "bonferroni")

bonf.result <- data.frame(
  Attribute = res$Attribute,
  Level = res$Level,
  Estimate = point.e,
  SE = sd.errors,
  p_raw = p.raw,
  p_bonf = p.bonf
)
bonf.result

plot.data <- bonf.result |> 
  mutate(
    signif = ifelse(p_bonf < 0.05, "significant", "not-significant")
  )

ggplot(plot.data, aes(x = Estimate, 
                      y = interaction(Attribute, Level, sep = " : "), 
                      color = signif)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = Estimate - 1.96 * SE,
                     xmax = Estimate + 1.96 * SE),
                 height = 0.2) +
  scale_color_manual(values = c("significant" = "red", "not-significant" = "gray40")) +
  labs(
    x = "AMCE",
    y = "Attribute: level",
    color = "Bonferroni",
    title = "AMCE（Bonferroni）"
  ) +
  theme_minimal(base_size = 16)


# Benjamini-Hochberg procedure
p.bh <- p.adjust(p.raw, method = "BH")
print(round(p.bh[p.bh <= 0.05], 3))

bonf.result$BH <- p.bh
bonf.result

plot.data <- bonf.result |>
  mutate(
    BH = p.bh,
    signif_BH = ifelse(BH < 0.05, "significant", "not-significant")
  )

ggplot(plot.data, aes(
  x = Estimate,
  y = interaction(Attribute, Level, sep = " : "),
  color = signif_BH
)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(
    xmin = Estimate - 1.96 * SE,
    xmax = Estimate + 1.96 * SE
  ), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("significant" = "red", "not-significant" = "gray40")) +
  labs(
    x = "AMCE",
    y = "Attribute: level",
    color = "BH(FDR)",
    title = "AMCE (Benjamini–Hochberg)"
  ) +
  theme_minimal(base_size = 16)



# Adaptive Shrinkage
result.summary <- summary(cjoint_pool)$amce

ash.norm <- ash(point.e, sd.errors, mixcompdist = "normal")

ash.point <- ash.norm$result$PosteriorMean
names(ash.point) <- paste(result.summary$feature, result.summary$level)

ash.CI <- ashci(ash.norm)
rownames(ash.CI) <- paste(result.summary$feature, result.summary$level)


ash.df <- data.frame(
  Attribute = result.summary$Attribute,
  Level = result.summary$Level,
  PosteriorMean = ash.norm$result$PosteriorMean,
  CI_low = ashci(ash.norm)[,1],
  CI_high = ashci(ash.norm)[,2]) |> 
  mutate(
    Label = paste(Attribute, Level, sep = " : "),
    Label = factor(Label, levels = rev(Label))
  )

ggplot(ash.df, aes(x = PosteriorMean, y = Label)) +
  geom_point(size = 3, color = "black") +
  geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height = 0.2, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  labs(
    x = "Posterior Mean (ASH)",
    y = "",
    title = "AMCE - Adaptive Shrinkage (ASH)"
  ) +
  theme_minimal(base_size = 14)

